#################################################################################
### Title: Are People Willing to Trade Away Democracy for Desirable Outcomes? ###
### Authors: Jonathan A. Chu, Scott Williamson, Eddy S. F. Yeung              ###
### Content: Subgroup analysis by country characteristics                     ###
### Date: September 20, 2025                                                  ###
#################################################################################

### Set-up ----
## Clean the working environment and set the working directory
rm(list = ls())
setwd("~/Desktop/democracy_tradeoff/replication") # change to your own working directory

## Load the required packages
library(tidyverse)
library(estimatr)
library(cregg)
library(expss)
library(cowplot)
library(grid)
library(gridExtra)
library(knitr)

## Read the cleaned and recoded datasets
df_US <- read.csv("df_US.csv")
df_JP <- read.csv("df_JP.csv")
df_EG <- read.csv("df_EG.csv")
df_IN <- read.csv("df_IN.csv")
df_IT <- read.csv("df_IT.csv")
df_TH <- read.csv("df_TH.csv")

## Merge the datasets
df_cj <- bind_rows(df_US, df_JP, df_EG, df_IN, df_IT, df_TH)

## Reorder the factors
# Leader selection
df_cj$`Leader Selection` <- 
  factor(df_cj$Leader.Selection,
         levels = c("Free and fair elections", "Unfair elections", "Unelected elites", "Hereditary succession", "Military coup"))

# Leader selection (binary)
df_cj$free_election <- ifelse(df_cj$Leader.Selection == "Free and fair elections", "Free and fair elections", "No free and fair elections")
df_cj$free_election <- 
  factor(df_cj$free_election,
         levels = c("Free and fair elections", "No free and fair elections"))

# Civil liberties
df_cj$`Civil Liberties` <- 
  factor(df_cj$Civil.Liberties,
         levels = c("Free", "Partially free", "Repressed"))

# Institutional checks
df_cj$`Institutional Checks` <- 
  factor(df_cj$Leader.Constraints,
         levels = c("Constrained", "Partially constrained", "Unconstrained"))

# National economy
df_cj$`National Economy` <- 
  factor(df_cj$National.Economy,
         levels = c("High income", "Middle income", "Low income"))

# Respondent wealth
df_cj$`Respondent Wealth` <- 
  factor(df_cj$Respondent.Wealth,
         levels = c("Wealthy", "Average", "Poor"))

# Public safety
df_cj$`Public Safety` <- 
  factor(df_cj$Public.Safety,
         levels = c("Very safe", "Somewhat safe", "Somewhat dangerous", "Very dangerous"))

# Corruption in politics
df_cj$`Corruption in Politics` <- 
  factor(df_cj$Corruption.in.Politics,
         levels = c("Rare", "Occasional", "Prevalent"))

# Health care
df_cj$`Health Care` <- 
  factor(df_cj$Health.Care,
         levels = c("Mostly accessible", "For the privileged"))

# Minority treatment
df_cj$`Minority Treatment` <- 
  factor(df_cj$Minority.Treatment,
         levels = c("Fairly treated", "Sometimes unfair", "Mostly unfair"))

# Respondent identity
df_cj$`Respondent Identity` <- 
  factor(df_cj$Respondent.Identity,
         levels = c("Majority", "Second largest", "Minority"))

# Country
df_cj$country <- 
  factor(df_cj$country, 
         levels = c("US", "JP", "EG", "IN", "IT", "TH"),
         labels = c("United States", "Japan", "Egypt", "India", "Italy", "Thailand"))

## Reverse the factors
df_cj$`Leader Selection` <- fct_rev(df_cj$`Leader Selection`)
df_cj$`Civil Liberties` <- fct_rev(df_cj$`Civil Liberties`)
df_cj$`Institutional Checks` <- fct_rev(df_cj$`Institutional Checks`)
df_cj$`National Economy` <- fct_rev(df_cj$`National Economy`)
df_cj$`Respondent Wealth` <- fct_rev(df_cj$`Respondent Wealth`)
df_cj$`Public Safety` <- fct_rev(df_cj$`Public Safety`)
df_cj$`Corruption in Politics` <- fct_rev(df_cj$`Corruption in Politics`)
df_cj$`Health Care` <- fct_rev(df_cj$`Health Care`)
df_cj$`Minority Treatment` <- fct_rev(df_cj$`Minority Treatment`)
df_cj$`Respondent Identity` <- fct_rev(df_cj$`Respondent Identity`)

## Bold feature labels in plots
bph <- c('bold', rep('plain', 5), 'bold', rep('plain', 3),
         'bold', rep('plain', 3), 'bold', rep('plain', 3), 
         'bold', rep('plain', 3), 'bold', rep('plain', 4),
         'bold', rep('plain', 3), 'bold', rep('plain', 2), 
         'bold', rep('plain', 3), 'bold', rep('plain', 3)) %>%
  rev() # reverse coding

## Create a function that indicates which estimates are statistically significant
## at the 5% level after using the BH procedure to adjust for multiple comparisons
sig1_fun <- function(data){
  data <- data %>% 
    mutate(sig = case_when(
      (p.adjust(p, method = "BH") < 0.05) == TRUE ~ 1, 
      (p.adjust(p, method = "BH") < 0.05) == FALSE ~ 0)) %>% 
    mutate(sig = factor(sig, levels = c(0, 1)))
  return(data)
}

### Asian countries vs. non-Asian countries (Figure S8) ----
df_cj <- df_cj %>% mutate(asian = case_when(
  country == "United States" | country == "Egypt" | country == "Italy" ~ 0,
  country == "Japan" | country == "India" | country == "Thailand" ~ 1,
))
df_cj$asian <- factor(df_cj$asian, 0:1, c("Non-Asian Countries", "Asian Countries"))

mms <- cj(data = df_cj,
          formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
            `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
            `Health Care` + `Minority Treatment` + `Respondent Identity`,
          id = ~id,
          estimate = "mm",
          by = ~asian)

p <- plot(mms, legend_pos = "none", size = 1.5) + 
  aes(shape = asian, color = asian) +
  xlab("Marginal Mean")
p <- p + scale_color_manual(values = c("darkred", "darkblue"))
p <- p + theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 9, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "gray50"),
    axis.title.y = element_text(size = 9, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )
p <- p + coord_cartesian(xlim = c(0.35, 0.65)) +
  geom_vline(xintercept = 0.5, color = "gray50", linetype = "dashed")
p + 
  annotate("text", x = 0.6299571, y = 40.25,
           label = "Non-Asian Countries",
           size = 2, color = "darkred", fontface = "italic", family = "Times") +
  annotate("text", x = 0.6417026, y = 41.75,
           label = "Asian Countries",
           size = 2, color = "darkblue", fontface = "italic", family = "Times")
ggsave("Figure S8.pdf", width = 6.5, height = 8)

## Present results in table form (Table S24)
table(df_cj$asian)
mms$N <- ifelse(mms$BY == "Asian Countries", 18426, 18474)
mms <- subset(mms, select = c("BY", "feature", "level", "estimate", "std.error", "lower", "upper", "N"))
kable(mms, digits = 3, "latex", row.names = F)

### Democracies vs. nondemocracies (Figure S9) ----
df_cj <- df_cj %>% mutate(democracy = case_when(
  country == "India" | country == "Thailand" | country == "Egypt" ~ 0,
  country == "United States" | country == "Japan" | country == "Italy" ~ 1
))
df_cj$democracy <- factor(df_cj$democracy, 0:1, c("Nondemocracies", "Democracies"))

mms <- cj(data = df_cj,
          formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
            `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
            `Health Care` + `Minority Treatment` + `Respondent Identity`,
          id = ~id,
          estimate = "mm",
          by = ~democracy)

p <- plot(mms, legend_pos = "none", size = 1.5) + 
  aes(shape = democracy, color = democracy) +
  xlab("Marginal Mean")
p <- p + scale_color_manual(values = c("darkred", "darkblue"))
p <- p + theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 9, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "gray50"),
    axis.title.y = element_text(size = 9, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )
p <- p + coord_cartesian(xlim = c(0.35, 0.65)) +
  geom_vline(xintercept = 0.5, color = "gray50", linetype = "dashed")
p + 
  annotate("text", x = 0.6268456, y = 40.25,
           label = "Nondemocracies",
           size = 2, color = "darkred", fontface = "italic", family = "Times") +
  annotate("text", x = 0.6447970, y = 41.75,
           label = "Democracies",
           size = 2, color = "darkblue", fontface = "italic", family = "Times")
ggsave("Figure S9.pdf", width = 6.5, height = 8)

## Present results in table form (Table S25)
table(df_cj$democracy)
mms$N <- ifelse(mms$BY == "Democracies", 18498, 18402)
mms <- subset(mms, select = c("BY", "feature", "level", "estimate", "std.error", "lower", "upper", "N"))
kable(mms, digits = 3, "latex", row.names = F)
